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We present detailed spectral and timing analysis of the hard x-ray transient 
IGR J16358-4726 using multi-satellite archival observations. A study of the 
source flux time history over 6 years, suggests that lower luminosity transient 
outbursts can be occuring in intervals of at most 1 year. Joint spectral fits of 
the higher luminosity outburst using simultaneous Chandra/ ACIS and INTE- 
^ GRAL/lSGRl data reveal a spectrum well described by an absorbed power law 

model with a high energy cut-off plus an Fe line. We detected the 1.6 hour pulsa- 
tions initially reported using Chandra/ ACIS also in the INTEGRAL/lSGKl light 
curve and in subsequent XMM-Newton observations. Using the INTEGRAL data 
^ ' we identified a spin up of 94 s (P = 1.6 x 10 -4 ), which strongly points to a neu- 

^ ; tron star nature for IGR J16358-4726. Assuming that the spin up is due to disc 

accretion, we estimate that the source magnetic field ranges between 10 13 ~ 10 15 



G, depending on its distance, possibly supporting a magnetar nature for IGR 
J16358-4726. 
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1. Introduction 

IGR J 16358-4726 belongs to a new class of transient hard X-ray sources discovered 
with ESA's INTErnational Gamma-Ray Astronomy Laboratory [INTEGRAL; Winkler et 
al. 2003). Kuulkers (2005) describes the spectral and temporal properties of ten of these 
sources discovered during regular scans of the Galactic plane with the Imager on-board 
INTEGRAL (IBIS; Ubertini et al. 2003) through the end of 2004. They all exhibit very 
high intrinsic photo-absorption (Nh > 10 23 cm -2 ) resulting in hard spectra with power law 
spectral indices 0.5 < V < 2.1. Four of these new transients exhibit long period coherent 
pulsations in their high energy (X- and 7— ray) fluxes in the range of 228 s to ~ 5900 s; 
so far they have been interpreted in the literature as High Mass X-ray Binaries (HMXBs), 
based on their spectral similarities to other X-ray pulsars in these systems (Kuulkers 2005; 
Lutovinov et al. 2005). In only one of these sources, however, an orbital period of 8.9 ±0.1 d 
has also been found, supporting an HMXB nature for the source (IGR J16320-4751; Corbet 
et al. 2005). 

Thus far, IGR J16358-4726 stands out as the sole outlier among these new INTEGRAL 
X-ray transients, having one of the hardest X-ray spectra and the longest period of all. 
The source was discovered with INTEGRAL by Revnivtsev et al. (2003) when it went into 
outburst on 2003 March 19. Our contemporaneous Chandra observations led to the discovery 
of large amplitude coherent intensity modulations in the source light curve with P = 5880 
s (Kouveliotou et al. 2003; Patel et al. 2004). Patel et al. (2004) showed that the source 
spectrum is extremely hard (r = 0.5) and highly absorbed (Nh = 1.3 x 10 23 cm -2 ), both 
during the outburst and in quiescence. Radio follow up observations taken ~ 25 days and 
~ 134 days after the outburst resulted only in flux upper limits of 0.9 mJy (4.8 GHz; Patel 
et al. 2004) and 7.5 mJy (0.61 GHz; Pandey et al. 2006), respectively. Wachter et al. 
(2006) report an accurate location of the source derived by off-axis PSF modeling of the 
initial Chandra observation. Their infrared observations, taken during the decaying tail 
of the source outburst, reveal a counterpart candidate for IGR J16358-4726 different than 
the one discussed by Kouveliotou et al. (2003) and D'Amico et al. (2006); their data allow 
both for a Low Mass X-ray Binary (LMXB) or an HMXB system scenario (based only on 
the brightness of the counterpart). The candidate suggested by (D'Amico et al. 2006) was 
observed during source quiescence; based on the lack of variability in the light curve of 
this putative counterpart, D'Amico et al. (2006) supported an HMXB nature of the source. 
Finally, based on the INTEGRAL survey observations, IGR J16358-4726 is thought to be 
located in the Norma arm of the Milky Way, at distance ranges of 5 — 6 or 12 — 13 kpc, 
depending on which crossing of the arm with the line of sight is chosen (Lutovinov et al. 
2005). 
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To date, there is no conclusive evidence for the nature of the IGR J 16358-4726 period. 
Patel et al. (2004) suggested that the source is either a transient accreting pulsar with an 
unusually long spin period, or a transient Low Mass X-ray Binary (LMXB) with a short orbit. 
Unraveling the nature of IGR J 16358-4726 is essential, as this source provides an excellent 
avenue towards understanding the new class of INTEGRAL objects. Here, we present a 
study of the source flux history and spectra using combined X- and 7— ray observations with 
ASCA, BeppoSAX, INTEGRAL, Chandra and XMM-Newton. We describe the INTEGRAL 
and XMM-Newton observations and data reduction in Section 2. In Section 3 we construct 
the long term flux history of IGR J16358-4726 spanning ~ 6 years and search for evidence 
of possible orbital modulations; further we perform a timing analysis using INTEGRAL 
data to identify changes in the observed period. In Section 4 we perform both a wide band 
joint spectral analysis to better constrain the integrated source spectrum and phase resolved 
spectroscopy with Chandra to search for spectral variations and features. A discussion of 
our results and possible models on the nature of the source is presented in Section 5. 

2. Observations and Data Reduction 
2.1. INTEGRAL 

A significant portion of the observing time allocated to the INTEGRAL Core Program 
is dedicated to the study of hard X- and 7-ray sources in the Galaxy. This time is mainly 
divided in Galactic Plane Scans (GPS) and Galactic Center Deep Exposures (GCDE). The 
observing strategy consists in pointings lasting 2.2 and 1.8 ks for GPS and GCDE, respec- 
tively, that are spread in various grid patterns to ensure the best use of the coded-mask 
instruments. As IGR J16358-4726 is located near the galactic plane, it has been regularly 
monitored by INTEGRAL (Lutovinov et al. 2005). Hereafter, we will reference the INTE- 
GRAL observations by revolution (Rev) number (1 total revolution is completed every ~ 3 
days) and by individual science window (sew) number (1 sew ~ 2 ks). IGR J16358-4726 
was discovered as a transient source in Rev 52 and faded beyond detection a few days later. 
Since then, no other detection of this source has been reported with INTEGRAL. 

We have searched the INTEGRAL public database to identify additional periods of IGR 
J16358-4726 activity; our search included all Guest Observer data obtained from the start of 
the mission until 2003 Oct 16.3 (MJD 52928.3), as well as Core Program data obtained until 
2004 Dec 2 (MJD 53341.1). Our data set comprises four epochs when the Galactic center 
was deeply observed resulting in a total exposure time of 4.3 Ms. The average off-axis angle 
for all IGR J16358-4726 observations is larger than 10°. We have, therefore, only used data 
taken with the large FOV coded- mask instrument IBIS/ISGRI (20 keV-1 MeV, Ubertini et 
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al. 2003; Lebrun et al. 2003) on board INTEGRAL in this investigation. Table 1 presents the 
details of the source visibility periods, and, for each of them, the fraction of the time when 
the source was effectively observed in the IBIS partially-coded field of view (PCFOV) of 
29° x 29°. Revolutions with net source exposure < 1 ks were discarded since these typically 
correspond to pointings during which the source crossed the PCFOV far from the center. 
These omitted revolutions are: 56, 61, 63, 164, 169, 179, 183, 241, 243 and 244. 

The data from all pointings with the source within the ISGRI PCFOV were reduced with 
the Offline Scientific Analysis (OSA) software 1 (v4.2 and v5.0) provided by the INTEGRAL 
Science Data Center (Courvoisier et al. 2003). This resulted in ~ 2000 science windows 
distributed between revolutions 30 and 234. For each pointing we extracted sky images and 
we then built mosaic images with longer exposures. As the source is a transient, we extracted 
and inspected light curves from each individual pointing as well as from the composite images. 
During the times the source was detected we extracted light curves on shorter time scales 
of ~ 50 s. We also generated individual spectra for each pointing and then summed them 
to obtain one average spectrum per revolution where the source was detected. The same 
technique was used to generate phase-resolved spectra based on the pulse period derived 
with the ISGRI data. All the pointings were initially analyzed with OSA 4.2, however, the 
revolutions in which the source was detected were subsequently re-analyzed with OSA 5.0 
(see details in Table 2). Lightcurves were corrected to the solar system barycenter using the 
OSA 5.0 barycent tool. 

It was practically impossible to observe the source with the X-ray instrument JEM-X, 
which has a FOV with a diameter of 13.2°; moreover, it was very difficult to analyze sources 
lying more than 3° away from the center of the JEM-X FOV. In the single pointing where the 
source was located ~ 2.5° from the center (revolution 54), we do not detect it with a 3cr upper 
limit of 0.005 count s -1 . Finally, data from the third instrument on board INTEGRAL, the 
Gamma Ray Imaging Spectrometer (SPI), were also not analyzed because the instrument 
angular resolution of 2° is too large for the crowded region where IGR J16358-4726 is located 
(e.g., 4U 1630-47 is located only 19' away). 

2.2. XMM- Newton 

We observed the location of IGR J16358-4726 with XMM-Newton starting on 2004 
February 15.6 UT and again on 2004 September 4.8 UT for ~ 32 ks in each pointing. The 
field of IGR J 16358-4726 was observed again by XMM-Newton during a scheduled observation 
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of SGR 1627-41 on 2004 September 23.3 UT. During the latter observation, the EPIC-PN 
was operating in Small Window Mode thus no PN data of IGR J 16358-4726 were collected. 
IGR J16358-4726 was detected 9'. 8 off-axis with the MOS detectors during this observation 

The observation data files were processed using the XMM-Newton Standard Analysis 
System (SAS v6.5.0) tools emchain and epchainaad HEASOFT/FTOOLS (v6.0.2). We first 
identified and removed times of extremely high background rates in the PN and both MOS 
detectors (> 15 counts s _1 across the entire detector) and then repeatedly applied a 3 o 
clipping criterion to the 0.5 — 10.0 keV count rates until all the remaining events were within 
±3 er of the mean rate. The remaining useful exposure times are given in Table 1. 

For the first two observations, when the source was detectable and imaged on-axis, we 
extracted source spectra from each detector using an extraction radius of r = 18"; we used 
r = 30" in the third observation when IGR J 16358-4726 was imaged off-axis. We then 
built instrument response files using the standard SAS tools arfgen and rmfgen. Background 
spectra were collected using offset regions close to the source location for each observation 
with a total area 15 times the source region. Finally, we grouped spectral bins to obtain at 
least 15 counts per bin before background subtraction in each spectrum. 

We extracted events centered on the Chandra source location from the first XMM- 
Newton PN observation; the source was not detected. Using WebPimms and assuming the 
Chandra spectral model (r = 0.8, log N u = 23.3; Patel et al. (2004)) we derived a So- 
upper limit on the 1 — 10 keV unabsorbed flux of 5.4 x 10~ 14 ergs cm" 2 s -1 . The limits 
derived from the MOS detectors are consistent with but not more constraining than the PN 
measurement. During the second XMM-Newton observation we clearly detected a source at 
the location of IGR J16358-4726 with an average total PN count rate of (1.78 ±0.09) x 10~ 2 
count s" 1 (1 — 10 keV). We estimated a PN background rate of (0.31±0.01) x 10~ 2 count s -1 . 
For the MOS1 and MOS2 detectors we found total rates of (5.5 ± 0.5) x 10~ 3 count s _1 
and (6.0 ± 0.5) x 10~ 3 count s -1 and calculated background rates of (1.06 ± 0.03) x 10 -3 
count s _1 and (1.53 ± 0.04) x 10 -3 count s -1 respectively. Source and background spectra 
were also extracted from the off-axis data from the final XMM-Newton observation. The MOS 
detectors clearly detect the source with total and background rates of (6.9 ± 0.4) x 10~ 3 
count s" 1 and (2.89 ±0.07) x 10" 3 count s _1 for MOS1 and (6.1 ±0.4) x 10" 3 count s" 1 and 
(3.74 ± 0.08) x 10~ 3 count s _1 for MOS2. However, the poor quality of the MOS spectra in 
this observation did not allow for useful spectral analysis. 
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3. Results 

3.1. Source position 

IGR J 16358-4726 was brightest during revolution 52 with an average count rate of 
5.8±0.1 count s _1 and a detection significance of 44.7a (Figure 1; 20 — 60 keV). We used a 
mosaic of that epoch to extract the INTEGRAL best source position at RA = 16 /l 35 m 52.5 s 
and Dec. = — 47°25'16" (J2000) with an uncertainty of 20" that includes a 10" systematic 
error due to instrument misalignment. An additional systematic error (proportional to the 
source significance) of 44" is also present from the image reconstruction (Gros et al. 2003). 
This position is consistent with the location initially reported by Kouveliotou et al. (2003b). 
A more accurate location using the second XMM-Newton observation and the initial Chandra 
observations is presented in a companion paper by Wachter et al. (2006) together with a 
discussion on the results of our search for an infrared counterpart for IGR J16358-4726. 



3.2. Source X-ray Light Curve 

Figure 2 exhibits the ISGRI light curve of IGR J16358-4726 during the 2003 March 
outburst of the source (Revolutions 52 through 55, spanning a total of 9 days). The gaps 
in the Figure correspond to Revs 53 and 56, when the source was outside the IBIS PCFOV. 
During the source maximum (Rev 52) the average ISGRI count rate was 5.8±0.1 count s _1 
corresponding to 36 mCrab. The initial peak lasted less than 2 days; the source flux fell 
to 3.2 count s -1 (20 mCrab) during the following 8 days and beyond detection level < 1 
count s" 1 (< 7 mCrab) after Rev 57. 

The long term flux history of IGR J16358-4726 (unabsorbed 2 — 10 keV) spanning 
~6 years was constructed by using all available pointed observations of the source and is 
presented in Figure 3. The archival source flux measurements with ASCA and BeppoSAX 
prior to the 2003 outburst were originally presented in Patel et al. (2004). The 2 — 10 keV flux 
was estimated for each X-ray instrument by assuming a spectrum similar to the one measured 
with Chandra in 2003. The INTEGRAL flux was extrapolated using the best fit spectral 
model to the ISGRI data during the outburst (Rev 52). Using the extrapolated INTEGRAL 
discovery flux and the two Chandra flux measurements, we find that the outburst lightcurve 
decays as an exponential with an e-folding time of 9.2 ± 0.6 days. The RXTE observation 
of IGR J16358-4726 is complicated by significant X-ray emission from the Galactic ridge. 
Using spectral fits obtained after modeling and subtracting the Galactic ridge component 
(Revnivtsev et al. 2003), we estimate an unabsorbed 2 — 10 keV flux of 2.9 x 10~ 10 erg cm _2 s _1 
on MJD 52724.4. This flux is nearly a factor of 2 larger than the flux we derive from our 
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CXO measurement on MJD 52722.3. Such a large change in the source flux is diffcult to 
reconcile with the constant average rate observed during the 7.1 hour CXO observation. 
This data point has been omitted from Figure 3 for clarity. The corresponding 2 — 10 keV 
luminosities in Figure 3 were then determined by assuming an average source distance of 7 
kpc. 

The detections and non-detections of IGR J16358-4726 shown in Figure 3 motivated 
a search for orbital periods based on the assumption that this system exhibits a behavior 
similar to the largest class of known transient accreting pulsars, those with Be (or Oe) star 
companions. These systems often have series of outbursts separated by their orbital period, 
with the outbursts occurring near periastron passage in an eccentric orbit (Bildsten et al. 
1997). For each one of a finely spaced grid of trial orbital periods in the range of 10 — 2000 
days, an attempt was made to divide the orbit into an active phase region (containing all the 
detections), and a quiescent phase region (with all the non-detections). Each time such a 
phase segregation was possible, it was considered as a candidate period. Figure 4 shows the 
duration of the longest possible quiescent orbital phase region for each candidate period. Note 
that no candidate periods longer than ~ 400 days were found. Thus, although the search did 
not identify any unique orbital period, it did demonstrate that a, yet undiscovered, periodic 
outbursting could be occurring in intervals of at most 1 year. 



3.3. Timing analysis 

We searched for all detections of IGR J16358-4726 in the public, two-year long IBIS/ISGRI 
database (see also Section 2) in two energy ranges: 20 — 60 keV and 60 — 150 keV. We then 
built mosaic images for each revolution and extracted light curves at each energy band where 
the net exposures (at the source location) were higher than 5 ks, for better statistics. At 
no time was the source detected above 60 keV. Whenever the source was not detected in 
one revolution (typical 3a flux upper limits of ~ 0.3 count s" 1 ), we grouped adjacent revo- 
lutions to increase sensitivity. The final net exposure times for the four resulting detection 
revolutions and their dates are given in Table 2. 

Initially, we generated Lomb-Scargle periodograms (Scargle 1982) from the ISGRI light 
curves with the higher count rate detections (Revs 52, 54, and 55; Figure 2) and clearly 
detected significant modulations at ~ 1.65 h (as well as its first harmonic). Further, we 
divided the light curve into 2 sections (Rev 52 and Rev 54+55) to search for changes in 
modulation frequency. Although we detected significant pulsations in the initial portion of 
the light curve (until ~7x 10 4 s after detection) with a period of 5970 ± 15 s, in the latter 
portion of the light curve we found moderate evidence for pulsations consistent with the 



- 8- 



previously measured frequency. In both orbit segments we detected significant power in the 
~ 3000 s harmonic. However, due to the poor statistics during the second ISGRI interval, 
we were unable to conclusively determine pulse period changes exclusively within these data 
using the Lomb-Scargle method. Finally, we proceeded into folding the first ~ 10 5 s of the 
ISGRI light curve on the best determined INTEGRAL period (Figure 5) and compared it 
with the Chandra X-ray light curve acquired on 2003 March 24 folded on the same ephemeris. 
There appears to be evidence of a phase shift b etween these two profiles suggesting evolution 
of the pulse profile or period over a few days; however, we caution that since the profiles are 
in different energy bands this shift could be due to spectral evolution. 

In the following, we describe a more detailed method for the determination of the pulse 
ephemeris. We again determine the pulse periods separately from the data in Rev 52 and 
Revs 55+54 combined using a method that takes into consideration the properties of the data 
quality and the observation windows. An important aspect of this analysis was the need to 
correctly account for the data gaps and the large variations in rate errors, which occurred in 
both of these intervals. These gaps and variations were due to changes in pointing directions 
during the observations, causing missing data when the source was outside the coded field 
of view, and increasing the errors as the source moved outward within the PCFOV. For this 
purpose only exposure-corrected count rates in the 20 — 45 keV energy band with 50 s time 
resolution were analyzed. 

We first determined the pulse periods using a technique based on y 2 fits of the rates. To 
account for the gradual decrease in flux during the observations, the rates in each interval 
were fitted to separate linear trends, and the rates were then divided by the fit. The trend 
amounted to 1.0% in Rev 52, and 31% over the course of Revs 54+55. For a grid of trial 
pulse periods ranging from 5000 s to 7000 s, the normalized rates were then fitted to a 
model consisting of a repeated pulse profile, parameterized by a Fourier expansion limited 
to constant, fundamental and first harmonic terms. Most of the power of the pulse was 
found to be in the fundamental and first harmonic, so we did not include higher harmonics. 
The improvement in this fit (Ax 2 ) over a fit of the normalized rates to a constant is shown 
versus the period in Figure 6 for both time intervals. This A% 2 is equivalent to the Z\ 
statistic (Buccheri et al. 1983) in the case of uniformly spaced data with uniform errors and 
periods much smaller than the duration of the observation. The best fitted periods were 
P = 5964.7+ 12.1 s (Rev 52) and P = 5870.9 + 3.4 s (Rev 54+55), where the errors include 
only counting statistics. 

To check these results and to quantify any noise in excess of counting statistics, each 
analysis interval was divided into sub-intervals and the pulse profile for each sub-interval was 
obtained by fitting the rates with a second order Fourier expansion pulse profile model, with 
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the period fixed to that determined above for that interval. These were then cross-correlated 
with the mean profile for the full interval to obtain pulse phase measurements. The deviation 
of these phases from those expected from the estimated period is shown in Figure 7. Using 
linear fits to the pulse phases, we obtained pulse period estimates of P = 5964.9 ± 12.2 s for 
Rev 52, and P = 5871. 3±3. 3 s for Revs 54+55 combined, consistent with the previous results. 
The x 2 °f these fits is 4.72 and 3.47, respectively (both with three degrees of freedom). These 
results are summarized in Table 3 with the errors increased by the square root of the reduced 
chi-squares to account for any noise in excess of counting statistics. 

In Figure 8 we show the evolution of the IGR J16358-4726 pulse periods. The two 
Chandra measurements agree with the second ISGRI measurement (Rev 54+55), but the 
initial ISGRI period (Rev 52) is significantly higher than the others. The 94 s difference 
between the first and second ISGRI periods is significant at the 6cr level. If this difference is 
due to a gradual change in the spin period of a neutron star, the mean spin period derivative 
between the measurements is ( — 1.65 + 0.28) x 10 -4 . The last Chandra measurement indicates 
that this putative spin-up trend may have ended. Unfortunately, the Rev 54+55 data only 
weakly constrain period changes; fitting the pulse phases in Revs 54+55 with a quadratic, 
results in the estimate of P = (1.9 ±1.1) x 10~ 4 . Over the 18.8 pulse cycles which occur 
during Rev 52, the 94 s difference in period amounts to only 1760 s or 0.30 pulse cycles of 
phase difference. It could be possible that the pulse profile, which our analyses assume is 
constant, has changed sufficiently to cause such an apparent phase shift. Figure 9 shows the 
20 — 45 keV ISGRI rates epoch-folded with a 5965 s period, for the five sub-intervals of Rev 
52 used in Figure 7. Also shown on the profile is the position of the phase reference that 
would have resulted if the profiles were folded with a period of 5871.3 s. Examining these 
alternate phase alignments of pulse profiles, we conclude that this period difference is most 
likely not due to profile changes. 

We have also determined a pulse period using the XMM-Newton EPIC-PN observation 
from 2004 September; in all other data sets the source was too faint to detect pulsations. 
We correct the photon times of arrival in the events list to the solar system barycenter and 
extract all events between 0.5 — 10.0 keV from a slightly smaller extraction region (12") as 
was used in the spectral analysis (see also Section 4). We performed a Z\ test on periods 
ranging from 4000 — 7000 sec and clearly detected the pulsation in the EPIC-PN light curve 
with a period of P = 5858 ± 74 s. This period is consistent with those measured using 
Chandra and INTEGRAL in 2003. We present the folded PN light curve in Figure 10. 
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4. Spectral analysis 

4.1. INTEGRAL /ISGKL + Chandra/ ACIS-S Spectral Analysis 

The phase averaged spectral fit results from the Chandra data alone have already been 
presented in Patel et al. (2004). We have, therefore, initially fitted a power law model 
to the INTEGRAL/ISGKL 20 - 200 keV spectrum, which resulted in an unacceptable fit 
l v = 79/9). A thermal blackbody (BB) model resulted in a better fit (x 2 /v = 13.7/9) 
with kTsB — 6.2 ± 0.1 keV; however, the resulting BB source radius was only 0.15 km 
(assuming a source distance of 7 kpc). Subsequently, we jointly fitted the phase averaged 
spectra measured with ISGRI and the initial Chandra observation, to better characterize the 
IGR J 16358-4726 spectrum in outburst. Since these data were collected at slightly different 
times, we allowed the model normalizations between data sets to be free and linked all other 
spectral parameters. The data were well fit by a Comptonization model (Titarchuk 1994) 
with the best fit parameters given in Table 4. For comparions to power-law models typically 
used in the lower x-ray bands, we have also modeled the spectrum using an absorbed power- 
law model with a high energy exponential cut-off (Table 4). 

On 2003 March 24, we were fortunate to have simultaneous INTEGRAL and Chandra 
observations providing us with the only broadband high energy spectral coverage for IGR 
J16358-4726 to date. The INTEGRAL data available during the time of the Chandra obser- 
vations consist of 4 science windows from which we were able to extract ~ 4600 s of useful 
exposure time. The light curves from these observations are shown in Figure 11; the overall 
shape of the ISGRI light curve is consistent with the shape of the Chandra light curve during 
the same time interval suggesting no evolution in the hard x-ray ephemeris. We performed 
joint phase-averaged spectroscopy using these simultaneous observations covering approxi- 
mately one entire pulse (Figure 11). The data are well represented by a Comptonization 
model with a Gaussian line feature at 6.4 keV. The best fit parameters are given in Table 4 
and the resulting spectrum is shown in Figure 12. 

4.2. XMM-Newton/EPlC-PN Spectral Analysis 

We employed standard x 2 methods in obtaining the best fit to the XMM-Newton spec- 
tral data and analyzed the spectra using the XSPEC (vl 1.3.21) spectral-fitting package (Ar- 
naud 1996). When fitting the MOS and PN spectra simultaneously we linked the model 
parameters. However, since the relative flux calibration of these instruments is somewhat 
uncertain, we allowed the individual normalizations to be free to obtain their best fit value. 
We fitted the MOS and PN (0.5 — 10.0 keV) spectra jointly to a power-law and again to a 
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thermal BB model. We accounted for absorption due to intervening gas and dust from the 
Galaxy using the Tuebingen-Boulder ISM absorption model (XSPEC model tbabs) with the 
relative abundances and cross-sections from Wilms et al. (2000) and with He cross section 
from Yan et al. (1998). An acceptable fit (x 2 l v — 30.3/38) resulted when the data were 
modeled with the power law with a best fit spectral index T = I.II03 an d column density 
Ah = 18l 3 x cm - 2 - The unabsorbed (1 — 10 keV) flux measured from this model was 
7.1 x 10~ 13 ergs s" 1 cm" 2 . A fit using a BB model resulted in a nearly identical quality fit 
(x>/v = 30.4/38) with kT BB = 2.3+°;^ keV and A H = (14 ± 3) x 10 22 cm' 2 . The unabsorbed 
(1 — 10 keV) flux measured for this model was 4.9 x 10" 13 ergs s™ 1 cm" 2 . 

To test the goodness of the x 2 fitting, we also fitted the unbinned data using the C- 
statistic Cash (1979) and again using the grouped data and the x 2 fitting with the weighting 
method described by Churazov et al. (1996) as implemented in XSPEC, and found consistent 
results. 



4-2.1. Phase Resolved Spectroscopy 

We have extracted 10 individual phase resolved spectra from the Chandra ACIS obser- 
vation of 2003 March. We simultaneously fit the spectra using an absorbed power-law plus 
Gaussian line model. The model consists of 6 parameters for each spectrum (the equivalent 
hydrogen column (A H ), spectral power-law index (r), power-law normalization (PLNORM), 
Gaussian line centroid energy (Eu ne ), width (crime), and flux (Fu ne )) for a total of 60 param- 
eters. For all fits we allow PLNORM for each spectral model to aquire its best fit value 
(i.e. PLNORM is free and unlinked). All the remaining parameters were forced to have the 
same value (i.e. were linked) for each spectrum and each parameter was free to acquire its 
best fit value. This fit results in a statically acceptable fit with x 2 l v — 562.6/542 and best 
fit model parameters of A H = 3.5 ± 0.1 x 10 23 cm" 2 , T = 0.59 ±0.09, E Hne = 6.38 ±0.01 keV, 
cr line = 54 ± 18 eV, and line flux of F Une = 2.6 ± 0.3 x 10~ 12 ergs cm" 2 s _1 . 

To search for variations with phase of the linked parameters, we first unlinked the Ah 
while continuing to hold the others linked. This resulted in significant improvements in the 
fit with Ax 2 = 29.1 with an F-test probability of F prob = 0.0011. We next unlinked T and 
found an additional improvement of A% 2 = 21.9. To test whether variations in T or A H (or 
both) were responsible for the majority of the fit improvement we repeat the same test this 
time first unlinking T. When unlinking only T we find a Ax 2 = 41.8 with F prob = 7.5 x 10" 6 . 
Unlinking Ah at this point negligibly improves the fit with an additional A% 2 = 9.5. There 
is clear evidence for change of spectral index with phase with the spectral index ranging 
from ~ 0.1 —0.9. Figure 13 shows the pulse profile and the evolution of the power-law index; 
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the spectrum is harder during the peak of the pulse and softens as the rate declines. We 
conclude that the evidence for varying column density with phase during this observation is 
marginal at best. 

We note that in Patel et al. (2004), we reported a significant decrease in the measured 
N a of (by ~ 40%) between this Chandra pointing and the one 27 days later. In subsequent 
spectral fits the power-law index remains unlinked. 

The emission line at 6.4 keV is most likely due to Fe fluorescent emission. We tested for 
variation in the 6.4 keV line parameters using a similar strategy. We performed a series of fits 
of the spectra after unlinking each line parameter while holding the remaining parameters 
linked and also performed a fit will all line parameters unlinked. In all cases the improvement 
in the spectral fits were not significant (Ax 2 < 8 for 9 additional parameters). We find no 
evidence for evolution of the line with pulse phase. 

Finally we searched the individual spectra for additional line features. The most sig- 
nificant line-like feature in a single spectrum occurred during phase 0.2-0.3 and appeared as 
a marginally significant emission line feature at 7.0 keV (A% 2 = 14.5 for 2 additional free 
parameters and width fixed at <Ju ne 2 = 50 eV). However, when we considered the number 
of spectral phase bins searched, we concluded that the line was not statistically significant 
with a probability of chance occurance of > 50%. 

5. Discussion 

There are two feasible origins for the pulsations observed in the flux of IGR J16358-4726: 
binary orbital modulation or neutron star spin modulation. Given the substantial difference 
of the period measured in Rev 52 of the INTEGRAL data from all later measurements, and 
the expected stability of an orbital period, we must either reject the Rev 52 measurement 
as unrepresentative of the intrinsic period, or reject the orbital modulation hypothesis. It is 
also possible that pulse shape evolution during Rev 52 could have resulted in this discrepant 
period, however this evolution would need to mimic a gradual shifting in phase by 30% of 
a cycle over the course of the observation. It is not clear how this could be explained in 
the context of an orbital modulation. Finally, attributing the 1.65 hour modulation to an 
orbital period would presumably imply that the object is an LMXB making the nature of 
the source difficult to reconcile with its unusually hard spectrum. 

Assuming the pulsations are spin induced, there are two possible explanations for the 
change in spin-period between the first observation and the following observations. Either 
the neutron star was spun-up by accretion induced torques, or the neutron star had a glitch 
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in spin period, analogous to those seen in rotation powered pulsars. 

We first consider a scenario where the observed modulation is the neutron star spin and 
the abrupt change in the frequency is attributed to the sudden tranfer of angular momentum 
between the superfluid interior and the pulsar crust (i.e., a "glitch"; Anderson & Itoh 1975; 
Ruderman 1976; Alpar 1977). We measure a change of frequency of Av = (2.7±0.4)/xHz and 
thus a fractional change in the frequency of Avjv = 1.6~ 2 . This is ~ 10 4 times larger than 
the values derived from large glitches reported from radio pulsars (e.g., Shemar & Lyne 1996; 
Krawczyk et al. 2003) and anomalous x-ray pulsars which have typical Avjv ~ 10~ 6 (Heyl & 
Hernquist 1999; Kaspi et al. 2000, 2003; Morii et al. 2005). To date, there has been only one 
substantiated report of a glitch observed from an accreting pulsar (KS 1947+300, Be system; 
Galloway et al. 2004). The fractional change in frequency measured in the KS 1947+300 
glitch was 3.7 x 10~ 5 , also significantly larger than those seen in radio or anomalous X-ray 
pulsars, but still much smaller than observed in IGR J16358-4726. We conclude, therefore, 
that unless the source underwent a unique and new type of glitch, never observed before, we 
can exclude this phenomenon as the cause of the neutron star spin variation. 

In the following scenario, we will now assume that the neutron star is fed by an accretion 
disk during the outburst observed by INTEGRAL. Then from the measured spin-up rate 
and flux, and an assumed distance, d, we can determine the pulsars magnetic moment. If 
instead the pulsar is wind fed, the spin-up at a given mass accretion rate is always less than 
the spin-up rate in the disk fed case, so that the magnetic moment calculated for the disk-fed 
case will be a lower limit for the wind fed case. The expected spin-up rate for disc accretion 
is (see e.g. Pringle & Rees 1972) 

v = (27iI)- 1 r ] (GMr A ) 1/2 M (1) 

were I is the neutron star's moment of inertia, M its mass, M the accretion rate onto the 
neutron star, and r A the Alfven radius for spherical accretion, 

r A = (2GM)- 1 V /7 M _2/7 (2) 

where \i is the neutron star's magnetic moment. The model dependent factor r] is 
expected to be near one in the "slow rotator" regime, where the rotation rate of the neutron 
star is much smaller than the rotation rate of the inner disk. For example, in the model of 
Ghosh & Lamb (1979, eqns. 10 and 11) this factor is given by rj = n(0)C 1/2 (= 1.39X0.52 1 / 2 = 
1.00) for a slow rotator. 

From the 94 s period change between MJD 52717.84 and 52724.39 seen in the INTE- 
GRAL data, we infer an average spin-up rate of v = 4.7 x 10~ 12 Hz s _1 (P = 1.6 x 10~ 4 ). 
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For the same time interval we estimate an average unabsorbed 2-100 keV flux of F — 
7 x 10 -10 erg cm _2 s _1 , using flux estimates from the joint INTEGRAL-Chandra fit of 
1.2 x 10~ 9 erg cm _2 s~ 1 during Rev 52, and 5.0 x 10~ 10 erg cm _2 s~ 1 at the time of the 
initial Chandra observation, combined with the trend of the IBIS/ISGRI 20-45 keV rates. 
With these values we find from equations 1 and 2 the magnetic moment 

H = 1.8 x 10 3 V 7/ W4X~ 3 4 7 V G cm 3 (3) 

where m 1A = M/1.4M©, R 6 = i?/10 6 cm, J 45 = J/10 45 gm cm 2 , and d 7 = d/7 kpc. 
Figure 14 shows the magnetic moment (and field) as a function of distance to the source 
assuming canonical neutron star mass and radius. 

As the mass accretion rate drops, the disk inner radius (or the magnetospheric radius) 
r o = C r 'A will expand to the co-rotation radius r co = (GM) 1 / 3 (2n / P)~ 2 / 3 , at which point 
the accretion will halt because material at r which is forced to co-rotate with the magnetic 
field is flung-off by the centrifugal force (Illarionov & Sunyaev 1975). The transition to this 
"propeller regime" should occur when the accretion flux of IGR J16358-4726 decreases to 

F mm = 1.6 x 10- V^Wi^e 7 4^~ 14 erg cnrV 1 (4) 

IGR J16358-4726 was detected in the XMM observations with a 1-10 keV flux of 7 x 
10~ 13 erg cm _2 s _1 . The spectrum was hard, and pulsations clearly present, both indicating 
that the neutron star was still accreting. This gives us a lower limit on the distance to the 
source of 

d > 6.7( e /0.25) 1 / 14 7 ? - 1 / 2 C -V4 m V6 jR -V2 / i/2 kpc (5) 
where e is the fraction of the total flux in the 1-10 keV band. 

The neutron star in IGR J16358-4726 was most likely born with a much shorter period 
than currently observed, and has been spun down to its current state. Davies & Pringle 
(1981) discussed this spin-evolution in the context of wind accreting systems, and showed 
that obtaining long period systems such as IGR J16358-4726 requires a combination of high 
magnetic moment and low mass capture rates from the stellar wind. The impact of accretion 
disk formation on this evolution is unclear. If we suppose the source is currently near an 
equilibrium period, with intervals of spin-up balanced by intervals of spin-down, than an 
interesting question is how long are the spin-down intervals? If the source falls back into the 
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sub-sonic propeller regime between outbursts the spin-down rate will be (Davies & Pringle 
1981) 

Pprop = 7 x 10- 4 i4 Q mZll^ s yr- 1 (6) 

where /i 30 = /x/10 30 G cm 3 . So the time needed to recover the 94 s period decrease seen 
in the during the outburst is 

Trecur = 94s / P prop = 3.5rf' 'm^Rp^ 'd) 2 yr (7) 

so that if the distance to the source is greater than 8.4 kpc, we are not likely to see 
another outburst from it in the next 30 years. The accretion disk solutions of Rappaport, 
Fregeau & Spruit (2004) show spin-down at low mass accretion rates which scale as equation 
6, but with a factor of six lower rate. If in fact the recurrance time is this long, we consider 
the fact that the transient was detected in outburst for the first time during the initial year 
of INTEGRAL Galactic plane scans as extremely fortuitous. Therefore, if the source is 
detected in a similar state again within the next 3.5 years, it most likely is a transiently 
accreting magnetar at ~ 7 kpc with a field strength of 3.4 x 10 14 G. 

SKP acknowledges support from NASA grants NNG04GH33G and DD3-4023X. MDS 
has been supported by the Italian Space Agency grant I/R/046/04. DE acknowledges the 
Israel Science Foundation, the US - Israel Binational Science Foundation, and the Arnow 
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Table I. Log of IGR J 16358 -4726 Observations 



Observatory 


Start 


End 


Revolutions 


Exposure 3, 




(MJD) 


(MJD) 






(days) 


INTEGRAL 


52650 


52743 


30- 


60 


22.953 


INTEGRAL 


52859 


52919 


100- 


119 


32.944 


INTEGRAL 


53060 


53084 


167- 


174 


10.164 


INTEGRAL 


53236 


53263 


226- 


234 


13.411 


Chandra 


52722.1691 


52722.4939 






0.297 


Chandra 


52750.1484 


52750.7093 






0.545 


XMM-Newton 


53051.6064 


53051.9865 






0.193/0.231/0.205 


XMM-Newton 


53252.7663 


53253.1255 






0.263/0.288/0.299 


XMM-Newton 


53270.6780 


53271.2782 






- /0.411/0.414 



Note. — All public and Core Programme ISGRI data when IGR J16358-4726 was within the ISGRI 
PCFOV between 2003 Jan 11 (Rev 30) and 2004 Sep 15 (Rev 234) are included. 

a The net useful exposure times. INTEGRAL times reflect ISGRI pointed observations where 
IGR J16358-4726 fell within the PCFOV. Times from XMM-Newton are for the EPIC PN/MOS1/MOS2 
detectors. 



Table 2. Log of INTEGRAL /ISGRI Observations with Detections of IGR J16358-4726 



Revolution 


Start 


End 


Exposure 


e a 


<?det h 




(MJD) 


(MJD) 


ks 


deg 




052 


52717.3972 


52718.7003 


77.1 


11.4 (6.1-18.2) 


44.7 


054 


52722.2812 


52724.6796 


89.6 


10.2 (2.5-18) 


28.6 


055 


52725.0607 


52726.4123 


90.4 


11.3 (5.7-14.3) 


16.3 


057 


52732.5051 


52733.6327 


70.5 


11.3 (5.7-11.3) 


6.1 



a Average off-axis angle. Numbers in parenthesis indicate the range of off-axis angle. Rev 53 and 56 were 
not considered since the source fell outside the PCFOV. 

b The detection significance is derived from the mosaiced images for each revolution (20-60 keV). 
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Table 3. Pulse Periods Determined From ISGRI Data 



Revolutions 


Epoch a 


Period 


Error 




(MJD) 


(s) 


(s) 


52 


52717.84 


5965 


15 


54-55 


52724.39 


5871.3 


3.5 



a Thc weighted mean of the observation times using the inverse square rate errors as weights. 
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Table 4. 


INTEGRAL k Chandra Joint Spectral Fit Results 




Parameter 


Value a 


Value 13 


Unit 




MODEL 1: Absorbed Comptonization Plus Gaussian 


A H 

a 

kTinjection 

kTpi asma 

T 

PC 

r 2-WkeV 
Tfd 

n ACIS/ISGRI 


27.9 ±0.9 
6.39 ±0.01 
0.07 ±0.02 
2.1 ±0.1 
7.9 ±0.6 
3.8 ±0.7 

2.7 ±0.03 

i 07+0.05 

1 - 01 -0.08 

0.45 
352.0/324 


36 ±2 
6.38 ±0.02 
0fi+ 03 

U.UD_ 06 

4+ 3 -3 

8.2 ±0.5 
9.5 ± 1.1 

3.3 ±0.7 

1 0+1.7 

L - ii -0A 
1 

140.5/120 


xlO 22 cm" 2 

keV 

keV 

keV 

keV 

x 10~ 12 ergs cm" 
x 10~ 10 ergs cm" 


" 2 s- 1 
2 s- 1 


MODEL 2: Absorbed Power Law w/ High Energy Cut-Off Plus Gaussian 


A H 
a 

r 

E cu t 

Efoid 

r 2-WkeV 

Racis/isgri 
X 2 /v 


33 ±1 

6.39 ±0.01 

0.07 ±0.02 

0.53 ±0.08 

19±l 

12.2 ±0.6 

2.9 ±0.03 
1 71+0.02 

1 i_o.09 

1.12 

361/324 


34 ±2 
6.39 ±0.02 

o.oetg-j 

0-39i°; 08 
23 ±3 
16 ±2.0 
3.2 ±0.7 

1 fi2 +0.04 

1 

140.5/120 


xlO 22 cm" 2 

keV 

keV 

keV 
keV 

x 10 _ 12 ergs cm" 
x 10~ 10 ergs cm" 


2 s- 1 
2 s- 1 



a Fit results derived using all useful ISGRI and ACIS data. Errors reflect the la uncertainty. 

b Fit results derived using ~ 4.6 ks of ISGRI and ACIS data collected concurrently. Errors reflect the la- 
uncertainty. 

c Unabsorbed flux derived from the ACIS spectra. 

d R is the ratio of the normalizations of the primary model component derived from the ACIS and ISGRI 
spectra. The model normalizations for both data sets are forced to be linked for fits to the data collected 
concurrently. 
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Fig. 1.— INTEGRAL IBIS/ISGRI (20-60 keV) image of the field containing 
IGR J16358— 4726 taken on 2003 March 19 (Rev. 52). The source detection was significant 
at the 44.7 a level in 77.1 ks of useful exposure time. Sources detected in this portion of the 
sky are labeled in the image and the locations of SGR 1627—41 (triangle) and 4U 1630—47 
(square) are also identified. The latter 2 sources are not detected above 20 keV during this 
revolution. For reference, the image pixel size is 4.9' on a side. 




Fig. 2.— ISGRI lightcurve of IGR J16358-4726 (20-45 keV) integrated in time bins of 
1000 s, starting on 2003 March 19 and spanning ~ 9 days. Revs 53 and 56 were not used in 
this plot because the source fell outside the IBIS PCFOV (see also text and Table 2). The 
dashed lines indicate the time interval of the intital 26 ks Chandra observation. 
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Fig. 3. — The unabsorbed X-ray flux history of IGR J16358— 4726; upper-limits at the 90% 
significance level are indicated for 3 BeppoSAX measurements and 1 XMM-Newton measure- 
ment. The fluxes (left side of Y-axis) derived from the INTEGRAL data are extrapolated 
into the 2 — 10 keV energy range using the best joint Chandra - INTEGRAL spectral fit for 
Rev 52. The luminosities (right side of Y-axis) were calculated assuming an average distance 
of 7 kpc for the source. 
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Fig. 4. — Maximum duration of the quiescent orbital interval for possible orbital periods 
of the IGR J16358— 4726 system. Only the periods are shown, where the orbit could be 
divided into an active phase region (containing all the detections), and a quiescent phase 
region (containing all the non-detections). 
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Fig. 5.— Top: INTEGRAL ISGRI pulse profile (20-45 keV) folded with a period P = 5965± 
15 s. Bottom: Chandra ACIS (2—10 keV) pulse profile folded with the same period. Both 
profiles are referenced to the same date (MJD=52700). The Chandra data were acquired 
~ 5 days after the ISGRI detection. 
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Fig. 6. — The decrease in \ 2 between a fit of the normalized rates to a constant and a second 
order Fourier expansion pulse profile is shown versus period (solid curve) for the ISGRI data 
in Rev 52 (top panel) and Revs 54+55 combined (bottom panel). Also shown (dotted curves) 
is the same statistic after the best fitted profile (at the peak A% 2 period) was subtracted 
from the normalized rates. This shows that most of the power away from the main peak is 
caused by the pulsations at the peak A% 2 period. 
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Fig. 7. — The phase offsets of the pulse profiles from sub-intervals of Rev 52 (top panel) and 
Revs 54+55 combined (bottom panel) relative to phase ephemerides with constant periods 
of 5964.7 s and 5870.9 s, respectively. 
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Fig. 8. — The pulse periods determined from the ISGRI data (•) along with the measure- 
ments based on Chandra observations (*). 
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Fig. 9. — Pulse profiles for five sub-intervals in Revolution 52 obtained by epoch- folding the 
20 — 45 keV rates at a period of 5965 s. The dashed lines show the location of phase zero 
which would have resulted if the rates were folded with the 5871.3 s period determined for 
Revolutions 54+55. 
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Fig. 10.— Folded XMM-Newton EPIC/PN pulse profile (0.5 - 10 keV) with period P = 
5858 ± 74 s measured in 2004 September and referenced to the same epoch as in Figure 5. 
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Fig. 11.— Simultaneous Chandra ACIS-S2 (2 - 10 keV) and INTEGRAL ISGRI (20 - 40 
keV) lightcurve of IGR J16358-4726. 
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Fig. 12.- Joint Chandra ACIS-S2 and INTEGRAL ISGRI vF v spectrum of 
IGR J16358— 4726 during outburst on 2003 March 24. A highly absorbed Comptonization 
model gives the best fit to the data. A Fe K a line feature clearly seen at 6.4 keV during the 
24 March observation had faded below the detection capability of Chandra during the second 
observation on 2003 April 21 (Patel et al. 2004). The resisduals to the best fit reported in 
Table 4 are shown in the bottom panel. 
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Fig. 13. — Phase resolved spectral fit results derived with Chandra ACIS-S2 (2.0—10.0 keV) 
data of IGR J 16358 -4726 in outburst on 2003 March 24. The top panel shows the 2-10 keV 
observed pulse profile (with 30 phase bins) and the bottom panel shows the spectral powerlaw 
index as a function of pulse phase (in 10 phase bins). There is evidence for evolution in the 
spectral index with the hardest spectra occurring near the end of the pulse rise and a clear 
spectral softening during the pulse decline. 
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Fig. 14.— The IGR J 16358-4726 2-100 keV luminosity, mass accretion rate, Alfven radius, 
magnetic moment, polar magnetic field strength, and recurrence time are plotted as a func- 
tion of distance. The y-axis units for each curve are described in the legend. The curves 
assume canonical neutron star mass, radius, and moment of inertia of m\ A — 1, Rq — 1, and 
I45 = 1. The vertical solid line reflects the lower limit on the distance based on our 2004 
XMM measurement. 



